Nathan Trouvain
Inria - Mnemosyne
R4 - 9 novembre 2021
Série temporelle de Mackey-Glass
L'équation de Mackey-Glass décrit l'évolution de différent phénomènes biologiques, comme la quantité relative de globules rouges dans le temps.
Cette équation est une équation différentielle retardée :
$$ \frac{dP(t)}{dt} = \frac{a P(t - \tau)}{1 + P(t - \tau)^n} - bP(t) $$où $a = 0.2$, $b = 0.1$, $n = 10$. Le retard $\tau$ est de $17$ pas de temps. $\tau$ a une forte influence sur le caractère chaotique de la série temporelle modelée. $17$ donne des résultats plutôt chaotiques.
from reservoirpy.datasets import mackey_glass
from reservoirpy.observables import nrmse, rsquare
timesteps = 2510
tau = 17
X = mackey_glass(timesteps, tau=tau)
# rescale between -1 and 1
X = 2 * (X - X.min()) / (X.max() - X.min()) - 1
plot_mackey_glass(X, 500, tau)
Prédire $P(t + 10)$ connaissant $P(t)$.
from reservoirpy.datasets import to_forecasting
x, y = to_forecasting(X, forecast=10)
X_train1, y_train1 = x[:2000], y[:2000]
X_test1, y_test1 = x[2000:], y[2000:]
plot_train_test(X_train1, y_train1, X_test1, y_test1)
units = 100
leak_rate = 0.3
spectral_radius = 1.25
input_scaling = 1.0
connectivity = 0.1
input_connectivity = 0.2
regularization = 1e-8
seed = 1234
from reservoirpy.nodes import Reservoir, Ridge
reservoir = Reservoir(units, input_scaling=input_scaling, sr=spectral_radius,
lr=leak_rate, rc_connectivity=connectivity,
input_connectivity=input_connectivity, seed=seed)
readout = Ridge(1, ridge=regularization)
esn = reservoir >> readout
y = esn(X[0]) # initialisation
reservoir.Win is not None, reservoir.W is not None, readout.Wout is not None
(True, True, True)
np.all(readout.Wout == 0.0)
True
L'apprentissage est offline ("hors-ligne") : il n'a lieu qu'une seule fois, sur l'ensemble des données d'entraînement.
esn = esn.fit(X_train1, y_train1)
2000it [00:00, 12340.04it/s]
plot_readout(readout)
y_pred1 = esn.run(X_test1)
100%|███████████████████████████████████████| 500/500 [00:00<00:00, 9522.81it/s]
plot_results(y_pred1, y_test1)
Coefficient de détermination $R^2$ et erreur quadratique normalisée :
rsquare(y_test1, y_pred1), nrmse(y_test1, y_pred1)
(0.9998810506313958, 0.0026034595249913566)
Passons d'un horizon de prédiction de 10 pas de temps à un horizon de 100 pas de temps
x, y = to_forecasting(X, forecast=100)
X_train2, y_train2 = x[:2000], y[:2000]
X_test2, y_test2 = x[2000:], y[2000:]
plot_train_test(X_train2, y_train2, X_test2, y_test2)
plot_results(y_pred2, y_test2, sample=400)
Determination coefficient $R^2$ and NRMSE:
rsquare(y_test2, y_pred2), nrmse(y_test2, y_pred2)
(0.94803604102033, 0.05593776748240906)
units = 500
leak_rate = 0.3 # - leaking rate
spectral_radius = 0.99 # - rayon spectral
input_scaling = 1.0 # - facteur de mise à l'échelle des entrées (input scaling)
connectivity = 0.1 # - densité des connexions du reservoir vers lui même
input_connectivity = 0.2 # et des entrées vers le reservoir
regularization = 1e-4 # - coefficient de régularisation (L2)
seed = 1234 # reproductibilité
esn = reset_esn()
x, y = to_forecasting(X, forecast=1)
X_train3, y_train3 = x[:2000], y[:2000]
X_test3, y_test3 = x[2000:], y[2000:]
esn = esn.fit(X_train3, y_train3)
2000it [00:00, 8100.12it/s]
seed_timesteps = 100
warming_inputs = X_test3[:seed_timesteps]
warming_out = esn.run(warming_inputs, reset=True) # échauffement
100%|███████████████████████████████████████| 100/100 [00:00<00:00, 4247.18it/s]
nb_generations = 400
X_gen = np.zeros((nb_generations, 1))
y = warming_out[-1]
for t in range(nb_generations): # génération
y = esn(y)
X_gen[t, :] = y
X_t = X_test3[seed_timesteps: nb_generations+seed_timesteps]
plot_generation(X_gen, X_t, nb_generations, warming_out=warming_out,
warming_inputs=warming_inputs, seed_timesteps=seed_timesteps)
Apprentissage se déroulant de manière incrémentale.
Utilisation de l'algorithme FORCE (Sussillo and Abott, 2009)
from reservoirpy.nodes import FORCE
reservoir = Reservoir(units, input_scaling=input_scaling, sr=spectral_radius,
lr=leak_rate, rc_connectivity=connectivity,
input_connectivity=input_connectivity, seed=seed)
readout = FORCE(1)
esn_online = reservoir >> readout
outputs_pre = np.zeros(X_train1.shape)
for t, (x, y) in enumerate(zip(X_train1, y_train1)): # pour chaque pas de temps de la série :
outputs_pre[t, :] = esn_online.train(x, y)
plot_results(outputs_pre, y_train1, sample=100)
plot_results(outputs_pre, y_train1, sample=500)
esn_online.train(X_train1, y_train1)
pred_online = esn_online.run(X_test1) # Wout est maintenant figée
100%|███████████████████████████████████████| 500/500 [00:00<00:00, 8039.78it/s]
plot_results(pred_online, y_test1, sample=500)
Determination coefficient $R^2$ and NRMSE:
rsquare(y_test1, pred_online), nrmse(y_test1, pred_online)
(0.9998104511137069, 0.0032864753423817064)
features = ['com_x', 'com_y', 'com_z', 'trunk_pitch', 'trunk_roll', 'left_x', 'left_y',
'right_x', 'right_y', 'left_ankle_pitch', 'left_ankle_roll', 'left_hip_pitch',
'left_hip_roll', 'left_hip_yaw', 'left_knee', 'right_ankle_pitch',
'right_ankle_roll', 'right_hip_pitch', 'right_hip_roll',
'right_hip_yaw', 'right_knee']
prediction = ['fallen']
force = ['force_orientation', 'force_magnitude']
plot_robot(Y, Y_train, F)
from reservoirpy.nodes import ESN
reservoir = Reservoir(300, lr=0.5, sr=0.99, input_bias=False)
readout = Ridge(1, ridge=1e-3)
esn = ESN(reservoir=reservoir, readout=readout, workers=-1) # version distribuée
esn = esn.fit(X_train, y_train)
100%|███████████████████████████████████████| 1444/1444 [00:51<00:00, 27.86it/s]
res = esn.run(X_test)
100%|█████████████████████████████████████████| 361/361 [00:14<00:00, 25.08it/s]
plot_robot_results(y_test, res)
print("RMSE moyenne :", f"{np.mean(scores):.4f}", "±", f"{np.std(scores):.5f}")
print("RMSE moyenne (avec seuil) :", f"{np.mean(filt_scores):.4f}", "±", f"{np.std(filt_scores):.5f}")
RMSE moyenne : 0.1693 ± 0.10356 RMSE moyenne (avec seuil) : 0.1445 ± 0.15254
Les données peuvent être téléchargées sur Zenodo : https://zenodo.org/record/4736597
im = plt.imread("./static/canary.png")
plt.figure(figsize=(5, 5)); plt.imshow(im); plt.axis('off'); plt.show()
display(audio)
Plusieurs motifs temporels répétitifs differents à décoder : les phrases.
im = plt.imread("./static/canary_outputs.png")
plt.figure(figsize=(15, 15)); plt.imshow(im); plt.axis('off'); plt.show()
esn = esn.fit(X_train, y_train)
100%|███████████████████████████████████████████| 90/90 [00:09<00:00, 9.22it/s]
outputs = esn.run(X_test)
100%|████████████████████████████████████████| 10/10 [00:00<00:00, 29289.83it/s]
scores # pour chaque chant testé
[0.9506604506604507, 0.9087872559095581, 0.9531759739013625, 0.9421525697165245, 0.9662828947368421, 0.9536231884057971, 0.8867638129933212, 0.9411598302687412, 0.9209350292196631, 0.9014840543100726]
print("Précision moyenne :", f"{np.mean(scores):.4f}", "±", f"{np.std(scores):.5f}")
Précision moyenne : 0.9325 ± 0.02503
units = 100 # - nombre de neurones dans le reservoir
leak_rate = 0.3 # - leaking rate
spectral_radius = 1.25 # - rayon spectral
input_scaling = 1.0 # - facteur de mise à l'échelle des entrées
connectivity = 0.1 # - densité des connexions du reservoir vers lui même
input_connectivity = 0.2 # et des entrées vers le reservoir
regularization = 1e-8 # - coefficient de régularisation (L2)
seed = 1234 # reproductibilité
Le rayon spectral est la valeur propre maximale de la matrice des poids du réservoir ($W$).
states = []
radii = [0.1, 1.25, 10.0]
for sr in radii:
reservoir = Reservoir(units, sr=sr, input_scaling=0.1, lr=leak_rate, rc_connectivity=connectivity,
input_connectivity=input_connectivity)
s = reservoir.run(X_test1[:500])
states.append(s)
units_nb = 20
plt.figure(figsize=(15, 8))
for i, s in enumerate(states):
plt.subplot(len(radii)*100+10+i+1)
plt.plot(s[:, :units_nb], alpha=0.6)
plt.ylabel(f"$sr={radii[i]}$")
plt.xlabel(f"Activations ({units_nb} neurons)")
plt.show()
$-$ rayon spectral $\rightarrow$ dynamiques stables
$+$ rayon spectral $\rightarrow$ dynamiques chaotiques
Rayon spectral et Echo State Property : rayon spectral $\rightarrow$ 1 (assure que les états internes ne sont pas affectés par l'initialisation).
Il s'agit d'un coefficient appliqué à $W_{in}$, venant changer l'échelle des données en entrée.
states = []
scalings = [0.01, 0.1, 1.]
for iss in scalings:
reservoir = Reservoir(units, sr=spectral_radius, input_scaling=iss, lr=leak_rate,
rc_connectivity=connectivity, input_connectivity=input_connectivity)
s = reservoir.run(X_test1[:500])
states.append(s)
units_nb = 20
plt.figure(figsize=(15, 8))
for i, s in enumerate(states):
plt.subplot(len(scalings)*100+10+i+1)
plt.plot(s[:, :units_nb], alpha=0.6)
plt.ylabel(f"$iss={scalings[i]}$")
plt.xlabel(f"Activations ({units_nb} neurons)")
plt.show()
Correlation moyenne des activités des neurones du reservoir avec les entrées :
for i, s in enumerate(states):
corr = correlation(states[i], X_test1[:500])
print(f"ISS : {scalings[i]}, correlation moyenne : {corr}")
ISS : 0.01, correlation moyenne : -0.2752749355697493 ISS : 0.1, correlation moyenne : -3.07462701783721 ISS : 1.0, correlation moyenne : -11.893465621922127
L'input scaling peut aussi être utilisé pour ajuster l'influence de chaque donnée en entrée.
avec $\alpha \in [0, 1]$ et:
$$ f(u, x) = \tanh(W_{in} \cdotp u + W \cdotp x) $$states = []
rates = [0.02, 0.2, 0.9]
for lr in rates:
reservoir = Reservoir(units, sr=spectral_radius, input_scaling=input_scaling, lr=lr,
rc_connectivity=connectivity, input_connectivity=input_connectivity)
s = reservoir.run(X_test1[:500])
states.append(s)
units_nb = 20
plt.figure(figsize=(15, 8))
for i, s in enumerate(states):
plt.subplot(len(rates)*100+10+i+1)
plt.plot(s[:, :units_nb] + 2*i)
plt.ylabel(f"$lr={rates[i]}$")
plt.xlabel(f"States ({units_nb} neurons)")
plt.show()
Le leaking rate contrôle la "mémoire" de l'ESN. Il peut être vu comme l'inverse de sa constante de temps